Gene expression in male and female stickleback from populations with convergent and divergent throat coloration

Abstract Understanding of genetic mechanisms underlying variation in sexual dichromatism remains limited, especially for carotenoid‐based colors. We addressed this knowledge gap in a gene expression study with threespine stickleback. We compared male and female throat tissues across five populations, including two in which female red coloration has evolved convergently. We found that the expression of individual genes, gene ontologies, and coexpression networks associated with red female color within a population differed between California and British Columbia populations, suggesting differences in underlying mechanisms. Comparing females from each of these populations to females from populations dominated by dull females, we again found extensive expression differences. For each population, genes and networks associated with female red color showed the same patterns for males only inconsistently. The functional roles of genes showing correlated expression with female color are unclear within populations, whereas genes highlighted through inter‐population comparisons include some previously suggested to function in carotenoid pathways. Among these, the most consistent patterns involved TTC39B (Tetratricopeptide Repeat Domain 39B), which is within a known red coloration QTL in stickleback and implicated in red coloration in other taxa.

males-although the literature on female ornament evolution has expanded in recent years (Hare & Simmons, 2019;Kraaijeveld, 2014).
There are two major classes of explanation for conspicuous female ornaments. The first is that ornaments are advantageous to females, through one or more of inter-sexual, intra-sexual, or social selection (Hernández et al., 2021). The second, which was essentially the default explanation for decades and is sometimes known as the genetic correlation or correlated response hypothesis, is that female ornaments evolve as a response to selection on males, owing to shared genetic underpinnings (Hare & Simmons, 2019;Kraaijeveld, 2014;Lande, 1980). Color patterns, the focus of the present work, are popular targets for studying the evolution of both ornaments and sexual dimorphism, as well as for the investigation of diverse aspects of natural selection (Kettlewell, 1955;Olendorf et al., 2006;reviewed in Cuthill et al., 2017;Orteau & Jiggins, 2020). The evolutionary genetics of color patterns were investigated extensively in the 20th and early 21st centuries using mendelian and quantitative genetic approaches (e.g., studies reviewed in McKinnon & Pierotti, 2010). With the spread of genomic methods, progress has accelerated in elucidating how specific color genes and pathways, some sexually dimorphic, evolve in a variety of systems. Carotenoid-based ornaments are of particular interest because they have been the subjects of widely cited studies of sexual selection (Hill, 1991, Milinski & Bakker, 1990 and the focus of considerable theorizing, especially with regard to honest signaling and good genes hypotheses (e.g., Koch & Hill, 2018). Nevertheless, understanding of carotenoid evolutionary genetics remains limited and this is an area of active research (Orteau & Jiggins, 2020;Toews et al., 2017).
Orange-red throat coloration has now been reported for female stickleback, as well as males, from at least three localities. These include widely separated drainages in which female color likely evolved independently from an ancestral state, in anadromous and/or marine populations, in which only males are red (McKinnon et al., 2000;Von Hippel, 1999;Yong et al., 2013Yong et al., , 2018. The genetics of female red throat color are of particular interest because there is so far no evidence that females with red throats are favored by sexual or social selection (Wright et al., 2015;Yong et al., 2013Yong et al., , 2015Yong et al., , 2018. The genetic bases for convergent adaptations in sticklebacks have often been shown to result from reuse of the same genes and sometimes the same alleles (Bassham & Catchen, 2018;Jones et al., 2012;Kitano et al., 2019;Roberts-Kingman et al., 2021; but see Fang et al., 2020), but whether convergent female ornaments result from the same or different genetic mechanisms is little studied in stickleback or in any other taxa (but see Yassin et al., 2016). While it is commonly assumed that genes and developmental mechanisms underlying similar traits in males and females are the same, this need not be the case. Indeed, van der Bijl and Mank (2021) report numerous examples from studies of mouse gene knockouts of cryptic sex differences in genetic architecture. Surprisingly, concordant regulatory changes may lead to discordant genetic effects in sexually monomorphic as well as dimorphic traits. In humans, genes with similar expression patterns in males and females may be regulated by different transcription factors and networks in each sex (Lopes-Ramos et al., 2020).
In QTL analyses of a California stickleback population with redthroated females, up to three genome regions mediated variation in female red throat coloration. The regions associated with male color were closely collocated and possibly the same. In addition, one region almost entirely overlapped a QTL for spine color, which is present in both sexes. The observation of extensively overlapping QTL for males and females is consistent with the hypothesis that female red throats arise largely as a correlated response to selection on males . Studies with other populations have found additional regions of the genome to be associated with variation in male color (Malek et al., 2012) and suggested a single-locus genetic architecture for red versus black male nuptial coloration (Hagen & Moodie, 1979). Some melanin-related genes have also been identified and characterized (Hart & Miller, 2017;Miller et al., 2007) in stickleback and the molecular genetic basis of variation in cryptic striping has been elucidated (Greenwood et al., 2011(Greenwood et al., , 2012. In addition, correlations (in some cases genetic) have been documented between male red color and other traits including female mating preference (Bakker, 1993), female body condition (Bakker et al., 1999), aggressive behavior (Bakker, 1994;Wright et al., 2016) and male vision (Brock et al., 2018).
Here we report the results of a study of gene expression in the skin of male and female threespine stickleback from populations in which males, both sexes, or neither sex may possess carotenoid-based throat coloration. Four populations are streamand freshwater-resident, three from British Columbia and one from California. A fifth population, from British Columbia, possesses an anadromous life history. Red female throat coloration has almost certainly evolved independently and convergently in the California and British Columbia stream-resident populations (Yong et al., 2013(Yong et al., , 2018. We address three main research questions. First, for British Columbia and California populations, is expression of the same or different genes and networks associated with convergent female throat color? Second, is female red throat coloration associated with expression of the same genes and networks as male red throat coloration, as predicted by the genetic correlation hypothesis? Third, for which genes is expression most strongly correlated with color, and are these genes known or candidate pigmentation genes, and/or associated with previously identified color QTL? 2 | ME THODS

| Overview of analyses and contrasts
In our main analyses, we used DESeq2 (Anders & Huber, 2010) to: identify genes whose expression was correlated with female red coloration within populations; test whether the expression of the same genes was correlated with female color in different populations; ask if the genes showing expression correlated with female color also showed differential expression in analyses including both males and females. Because some red-associated genes may differ in expression mainly between populations, we conducted complementary analyses comparing red females with dull females of different populations in which dull females predominate. To address similar questions at a systems level, we used WGCNA to ask if gene networks showing color-correlated differential expression, in each population possessing red females, were present in the other population; we conducted an analogous analysis of the two sexes. In order to identify candidate genes for further study, we asked if any of the genes consistently associated with red coloration had been linked with coloration in other species, and if any were associated with stickleback coloration QTL detected in one of the same populations . As a first step before the main analyses, we characterized the overall structure of the data and broad patterns of variation, by population and sex, with a principal component analysis of gene expression.

| Study populations
Fish were collected from populations in British Columbia, Canada, and California, USA (Figure 1) 2.3 | RNA extraction, library preparation, sequencing, and generation of read counts Red coloration was assessed following protocols described in Yong et al. (2013) using an Ocean Optics Maya2000 Pro spectrophotometer.
Prior to collection of spectrometry data, subject fish were euthanized by MS-222 overdose. Immediately following these measurements, throat and telencephalon tissue (these results to be presented in a different MS) were removed and stored in RNAlater. "Throat tissue" specifically refers to ventral surface tissue between the opercula and immediately ventral to the ceratohyal cartilage, excluding tissue immediately ventral to the sternohyoideus; this includes (potentially) red-colored skin tissue on the underside of the head, detached from underlying cartilage. Additionally, standard length was measured and each fish was dissected to confirm sex. The remaining tissue samples Following sequencing, standard quality control (adapter removal and trimming of low-quality bases) was conducted using TrimGalore on default settings (Krueger, 2015;Martin, 2011). Reads were mapped to the stickleback reference genome (Jones et al., 2012) using Tophat version 2.0.13 (Trapnell et al., 2009). Tables of read counts, which were used as input for most downstream analyses, were generated using HT-Seq version 0.6.0 (Anders et al., 2015). They referenced Ensembl 79.
Gene annotations were from Ensembl 100 unless otherwise indicated.

| Color measurement and analysis
Reflectance spectra were processed with pavo (Maia et al., 2013) in the R environment (R Core Team, 2017) using a perceptual model of threespine stickleback vision (Stuckert, Drury, et al., 2019) to generate a measure of chroma, r.vec, which can be thought of as color intensity, and two indices of wavelength, h.phi and h.theta. The r.vec variable, which in our analyses is a measure of red intensity and will henceforth be referred to as chroma, was the intended focus but all three variables proved to be strongly correlated (see Results). Some heteroscedasticity was present in chroma data owing largely to limited variance in the consistently low values of the Bonsall population (note that males in this site appear to be melanic based on our field observations; we have observed some red expression in males from this drainage held in the laboratory for an extended period, but not in males in the present study). This was not completely eliminated by transformations and significance values were high and robust, so analyses of raw data, conducted using JMP Pro 14.1.0, are presented for ease of interpretation. Reflectance data were not available for one Little Campbell Stream male, but that male's RNAseq results were retained in most analyses as other key data were present.

| Differential expression analyses
To characterize major patterns of variation in expression data, prin- fold change, which for an experiment is log 2 (treated/untreated); for continuous independent variables, as in some of our analyses, the reported log 2 fold change is per unit of change in the continuous variable.
Unless otherwise indicated, alpha values for DESeq2 analyses have been corrected for multiple testing using the method of Benjamini and Hochberg (1995) and are denoted as "padj." We accepted the DESeq2 default independent filtering of lowly expressed genes, which optimizes the number of adjusted p-values, resulting in a padj of "NA" for some genes. Uncorrected significance tests are included in supplemental tables as a complement to corrected tests. Provisional characterizations of genes not named in Ensembl are presented when highly significant or otherwise of interest. Because some populations do not have red females and only one population lacks red males, our sample design is not balanced. We therefore conducted our analyses through focused comparisons rather than through a comprehensive model.

Stream populations
We first tested whether expression of the same or different genes was associated with convergent female red throat coloration within

| Comparisons of red females and red males
Next, we tested whether genes whose expression correlated with red in females (1B above) were also associated with red across the sexes, as predicted by the genetic correlation hypothesis.
Correlations between chroma and gene expression were analyzed for males and females together, first with a univariate analysis and then with sex included as a covariate, and the results compared with those for genes identified in 1B. The analysis was then repeated with the addition of an interaction term for sex*red chroma, to test for different relationships between chroma and gene expression in males versus females. These within-population analyses were conducted using DESeq2 1.32.0. In addition, we assessed whether in inter-population contrasts, males and red females of the same populations differed in expression of the same genes when contrasted with females from populations lacking red in females (as in 3 below).

| Gene expression comparisons of red females with dull females from other populations
Because genes mediating variation in a trait within and between populations are not necessarily the same, we also conducted an "inter-population" analysis comparing red females (dull females were excluded to ensure a meaningful contrast) from each of LC Stream and Matadero against females from each of two freshwater-resident populations in which red females are less common (Salmon River) or absent (Bonsall Creek). That is, we compared red females from LC Stream with dull females from Salmon and with those from Bonsall, and the same for Matadero females. While these contrasts are expected to reveal genes that differ by population independent of color, genes associated with color will also be encompassed, especially genes that differ in expression between red females and both of the dull female populations.

| Functional enrichment testing
Sets of genes differentially expressed with red coloration in females of each red female population were tested for functional overrepresentation using GOrilla (November 25-30, 2020), which tests rank-ordered gene lists for functional enrichment toward the top of the list (Eden et al., 2007(Eden et al., , 2009. To run these analyses, stickleback Ensembl gene stable IDs were converted to those for zebrafish, which for some genes resulted in longer lists owing to zebrafish possessing multiple members of a gene family that all correspond to a single stickleback gene stable ID. To avoid artifacts we retained only the first (ordered by name) in the list of zebrafish genes corresponding to a stickleback gene; these results are presented. We tested the robustness of the results obtained by instead retaining only the last gene in the zebra fish list and rerunning the analysis. Using default GOrilla settings, padj of 0.05, and the "Process" ontology, the results were identical for the number of Gene Ontology terms shared between analyses with Matadero and with LC Stream females.

| Candidate gene annotation and analysis
Because some pigmentation pathways may not be thoroughly covered by gene ontology analyses (Baxter et al., 2019), we also compared differentially expressed genes with a list of candidate pigmentation genes associated with pigmentation/color phenotypes. First, we compiled stickleback genes orthologous to 650 genes in a manually curated list associated with integumentary pigmentation phenotypes in humans, mice, and zebrafish (Baxter et al., 2019). We used Ensembl (98; Sep. 29, 2019) to search for stickleback orthologs, retaining those which Ensemble assigned high orthology confidence (1 rather than 0). Of the 619 zebrafish genes in the list, we retained 526 stickleback orthologs; of 30 mouse genes lacking zebrafish homologs, we retained four; one human gene with neither a zebrafish nor mouse homologue lacked a stickleback homologue.
As a supplement to this extensively curated list of 530 genes, we used stickleback orthologs of genes in a color gene list assembled by Stuckert, Moore, et al. (2019), derived from a broader array of taxa but less extensively annotated and vetted. We first used Ensembl to search for stickleback orthologs corresponding to gene names. To be thorough, we then used the list to search for human and zebrafish genes, in turn using those to search for stickleback orthologs. We retained those with a stickleback gene name which Ensembl assigned high orthology confidence (1 rather than 0); we merged the three lists and removed duplicate genes. We augmented this set with genes associated with carotenoid pigmentation based on Toews et al. (2017), again using Ensembl to identify stickleback genes, a total of 10, corresponding to carotenoid-related genes BCO1, BCO2, and CYP2J19 (though the latter is not expected to have direct homologues in fish); we used multiple stickleback genes from the same subfamilies where one to one orthology could not readily be established. TTC39B has been found to be pigment related (Ahi, Lecaudey, Ziegelbecker, Steiner, Glabonjat, et al., 2020;Hooper et al., 2019;Salis et al., 2019) Table S1). We describe genes from this list as candidate pigmentation genes. In addition, genes were assigned to QTL Bayesian credible intervals from Yong et al. (2016) using the assembly update of Glazer et al. (2015).

| WGCNA analyses
We conducted network analyses of RNAseq data using the R package WGCNA 1.68. Weighted gene correlation network analysis (WGCNA) is a systems biology method that can be used to find clusters, called modules, of genes with highly correlated expression (Langfelder & Horvath, 2008;Langfelder et al., 2011). We analyzed networks based on expression data transformed using DESeq2's variance-stabilizing transformation (Anders & Huber, 2010), with all genes not expressed in at least 50% of samples removed and a minimum module size of 30. We tested the correlation of modules with red intensity in R, using the method of Benjamini and Hochberg (1995) to correct for multiple tests.
As with DE-Seq analyses, we used WGCNA to test for similarities and differences between Matadero and LC Stream populations. To do this, we built separate networks for each population, including males and females. Due to small sample sizes for LC Stream, we included the adjacent (and similar in overall gene expression) Salmon River population for network construction. The resulting modules were tested for associations with color (in both sexes and in females only). All modules, including color-associated modules, were then tested for preservation of network structure among populations. It is important to keep in mind that both males and females were included in these WGCNA analyses, when interpreting relationships between WGCNA modules and traits in a single sex, and that sample sizes were relatively small for this method.
We also tested for network preservation between sexes to identify whether genes associated with color in females were preserved in males. To do this required larger sample sizes, so we built another gene co-expression network for all of the British Columbia fish together. Using these modules, we tested for module preservation of color-associated modules between males and females.

| Major axes of variation in gene expression
RNA sequencing yielded a mean of 19.1 million reads per throat sample. On average ~83% mapped to the genome, and about 36% of mapped reads could be assigned to annotated genes. Reads that were not counted either mapped to multiple places in the genome, mapped outside of known features (i.e., genes), or could not be unambiguously assigned to genes.
In a PCA of throat expression profiles, PC1 and PC2 explained 39% and 12% of the variance, respectively (Figure 3). Two-way  Table S3).
A gene ontology analysis also supports the distinctiveness of color-associated genes in each of LC Stream and Matadero populations (Tables S4 and S5). For example, the top three (Process) terms for LC Stream were proteasomal ubiquitin-independent protein catabolic process; macromolecule catabolic process; cellular macromolecule catabolic process. In contrast, the top three terms for Matadero were rRNA metabolic process; rRNA processing; nucleic acid metabolic process. No terms directly associated with pigmentation were significant.
WGCNA analysis of LC Stream fish resulted in 28 co-expressed gene modules. The module eigengene, the first principal component, showed a significant correlation with red chroma for three (Table   S6). No modules correlated significantly with red chroma for just females (Table S7). When LC Stream modules were tested for preservation against the Matadero data set, Zsummary values ranged from 0.24 to 31.2, with three under two and therefore not preserved, 12 between two and 10 and therefore preserved but not strongly, and 13 over 10 and thus strongly preserved (Langfelder et al., 2011; Table S8). All modules exhibited high "quality," a measure of internal preservation of a module calculated using subsets of the original data. With three LC Stream modules essentially not detectable in the Matadero data, there are clearly differences between populations in the organization of their gene expression networks. However, all modules with a significant association with color (both sexes together) were preserved.
In the reciprocal analysis of Matadero, 28 modules were generated, of which three showed significant associations with red chroma overall (Table S9). One of these modules was significantly associated with red chroma for females specifically, orangered3 ( Figure 4, Table S10; module names are arbitrarily assigned by software). When tested against LC Stream data, Zsummary preservation values ranged from −.36 to 46.1, with four not preserved, three preserved but not strongly, and the remaining 21 strongly preserved.
Among the modules not preserved in this analysis was orangered3, the Matadero module most significantly associated with color in females and also significant overall. Thus, the differences between the populations in coexpression patterns were confirmed. The other two color-associated modules were strongly preserved. Once again all modules exhibited high quality (Table S11).

| Do genes showing red chroma-correlated expression in females show a similar pattern across sexes?
Within population DESeq2 analyses did not provide strong evidence for or against the hypothesis that the same genes mediate variation in red chroma within and between the sexes. When both males and females were included in an analysis of the relationship between red chroma and differential expression for the Matadero population, 275 genes were significant, but these included only 11 of the 28 genes significant in the analysis conducted exclusively with females (Table S12). When a term for sex was added to the analysis along with chroma, just one gene was significant for red chroma, fgl1, fibrinogen-like 1, whereas 375 genes differed in expression between the sexes. Finally, when the interaction between sex and chroma was also included, it was significant for only one gene, hpdl, 4-hydroxyphenylpyruvate dioxygenase-like (padj = .036).
For the Little Campbell Stream population, when both males and females were included in an analysis of the relationship between red chroma and differential expression, 10 genes were significant (Table   S13), though these included only three of the 83 genes significant in the female-only analysis. Once again, expression was significantly

F I G U R E 4 Relationship between red
chroma and module eigengene value (first principal component) for "orangered3" in WGCNA analysis of Matadero population, n=14. Females in red, males in dark grey correlated with chroma for fewer genes when a term for sex was added, just five (Table S14), in contrast to 636 that differed between the sexes. When the interaction between sex and chroma was added, it was significant for 13 genes (Table S15), though only two of these were also significant in the analysis confined to females.
Thus, there was limited evidence from either population that individual genes showing chroma-correlated expression in females were similarly correlated with red chroma in males; but there was also little statistically significant evidence that gene expression-chroma correlations differed between the sexes. Sample sizes likely limited the power of these analyses.
To test whether co-expressed modules of genes that correlated with chroma in one sex were present in the opposite sex, we conducted WGCNA analyses, focusing on British Columbia stickleback.
Analysis of females generated 32 modules, of which eight were significantly correlated with chroma (Table S16). When tested against male data, all modules were preserved and 26 were strongly preserved (Table S17). Analysis of males led to 25 modules, of which two were significantly associated with chroma (Table S18). When tested against female data, all male expression modules were preserved and 20 were strongly preserved (Table S19). Thus, the modules associated with color in each sex were present in the other.

| Red females from Little Campbell and Matadero differ from each other and from dull 34 females of other populations in overall and candidate gene expression
We evaluated differential expression for each of LC Stream and Matadero red females, relative to dull females from two other populations (Bonsall and Salmon, Table 1; genes identified here were significantly different from females of both dull populations). These patterns of differentiation reflect overall population differences in gene expression, with Matadero fish showing generally greater expression divergence, as expected from the PCA (Figure 2). They also reveal clear population differences in candidate pigmentation gene expression. Differentially expressed candidate genes were largely different for Matadero and LCS females and also from those revealed by within population analyses (Table 1). For Matadero (Table   S20), 34 differentially expressed candidate pigmentation genes included bco2b (Figure 5), previously implicated in studies of fish carotenoid pigmentation, and bco1 from the same gene family, as well as TTC39B (Figure 6). Just three candidate genes, TTC39B, csf1ra, and tmem138, were differentially expressed by LC Stream red females relative to dull females of other populations (Table S21) (Table S20). Of the 54 loci uniquely differentially expressed by LC Stream red females when compared with dull females of other populations, 30 showed significantly different expressions between LC Stream and Matadero red females (Table S21).
For females of the LC Stream population, there was no overlap between the genes showing a correlation between expression and red chroma within LC Stream and those showing differential expression between LC Stream red females and dull females of two other populations (Table 1; Table S21). For Matadero, nine genes from the 28 that were correlated with red within that population were also associated with red between populations. None were on the list of candidate pigmentation genes (Table 1; Table S20).
To gauge the sex specificity of the above results, we also contrasted males from LC Stream and Matadero against dull females from Bonsall and Salmon populations. Approximately half of the genes differentially expressed in red LC Stream females were also significantly differentially expressed in LC Stream males, relative to dull females. The same was true for 80% of genes for Matadero fish (Table 1, Tables S20 and S21). Matadero males differentially expressed 29 of the same candidate pigmentation genes as females, including some with roles in carotenoid pigmentation (Table S20), notably bco1, bco2b ( Figure 5), and TTC39B. For LC Stream, just two candidate pigmentation genes, TTC39B and tmem138 (Crawford et al., 2017) were differentially expressed by both males and red females (Table S21). Thus, TTC39B was one of the seven genes, and the only candidate pigmentation gene, differentially expressed by red females and males of both populations relative to dull population females. ENSGACG00000010436, noted above, was also differentially expressed by both sexes and populations ( Figure 6).

| DISCUSS ION
We Previous studies of the genetic basis of stickleback color pattern evolution have generally focused on melanic traits (e.g. Greenwood et al., 2012;Hart & Miller, 2017), quantitative genetic approaches (e.g., Bakker et al., 1999), or mendelian analyses based on organismal phenotype (e.g., Hagen & Moodie, 1979; but see Malek et al., 2012), and males have been studied almost exclusively (but see Yong et al., 2016). To the best of our knowledge, ours is the first study to investigate genomic patterns of expression in stickleback skin toward a better understanding of the genetics of red coloration, and one of the first to use any approach to investigate the contributions of gene expression to variation in female coloration and stickleback dichromatism. In other studies of stickleback evolving convergently or in parallel as they colonize new habitats, exactly the same mutations have sometimes been recruited from standing genetic variation. Extensive non-replicated divergence has also been observed, however, including in gene expression data (Bassham et al., 2018;Fang et al., 2020;Jones et al., 2012;Kitano et al., 2019;Roberts-Kingman et al., 2021).  (Colosimo et al., 2005;Jones et al., 2012). Similarly, in a study of gene expression divergence between stream and marine sticklebacks from different continents, more extreme divergence was associated with longer isolation . Other potential explanations for the Matadero population's distinctive expression patterns include evolution from a relatively divergent marine ancestor or different selection pressures.

|
Regarding the latter possibility, some differences in environmental selection are to be expected in light of Matadero's location more than 1200 km south of the British Columbia populations.
Somewhat unexpectedly, genes whose expression correlated with red within a population were mainly different from those identified by analogous between population analyses. In particular, clear differences in the expression of carotenoid-related genes were observed between populations, yet such genes were not associated with variation in red chroma within populations. This result can be explained if some expression differences are necessary for red color in females, but different genes mediate whether, or when, a potentially red individual actually expresses such coloration. Such a pattern may be consistent with the plasticity in male nuptial coloration that has been documented in some lake populations in response to fine-scale variation in the light environment (Brock et al., 2017). At the extreme, plasticity could even explain the differences in which genes were correlated (

| Variation in expression and coexpression patterns by sex
Our results did not clearly indicate whether the genes showing correlated expression with color in females were associated with color across the sexes, the pattern expected if female coloration evolved as a correlated response to selection on males (Kraaijeveld, 2014;Lande, 1980), rather than independently from males through possibly distinct mechanisms (Lopes-Ramos et al., 2020;Van der Bijl & Mank, 2021). Within populations, individual genes exhibiting correlated expression with female red coloration showed inconsistent patterns when males and females were analyzed together, yet significant differences between females and males for chroma-gene expression correlations were also rare. Comparisons between populations with or without red females tended to exhibit more overlap between the sexes in the genes associated with red coloration, but overlap was still far from complete. In network analyses, sets of co-expressed genes were largely shared between the sexes, including those that correlated with red chroma; but in this analysis and the previous, sample sizes limit our power to identify sex-specific patterns.
Nevertheless, in a few cases, these data provide evidence that specific candidate pigment genes likely influence red chroma in both sexes (discussed below), and patterns at the system level appear to be stronger.

| Candidate genes differentiate populations with red females and overlap with known QTL
TTC39Bwas the only candidate pigmentation gene significantly associated with color, in inter-population analyses, for both sexes of both populations possessing red females. This lipoprotein coding gene is found in a QTL identified by Yong et al. (2016), and thus may also have contributed to QTL results obtained with one of our study populations. Other studies of gene expression in fish suggest a role for TTC39B in orange-red coloration. In the cichlid, Tropheus duboisi (Ahi, Lecaudey, Ziegelbecker, Steiner, Glabonjat, et al., 2020), and the clownfish, Amphiprion ocellaris (Salis et al., 2019), TTC39B has been found to be upregulated in yellow or orange skin, relative to white.
In three phylogenetically matched pairs of cichlid species, TTC39B was consistently upregulated in the species with red, rather than yellow, skin color (Ahi, Lecaudey, Ziegelbecker, Steiner, Goessler, et al., 2020). Beyond fish, correlations between color and TTC39B expression have also been observed in poison frogs (Stuckert et al., 2021 Additional genes, previously shown or suggested to affect color patterns, were significantly associated with red female color for just a single population. In particular, bco2b was highly significant in the interpopulation analyses for red Matadero females and males. In general, it was expressed at a higher level in the dull fish, as expected for this gene, which breaks down carotenoid pigments. Reduced bco2b expression has been linked with persistence and accumulation of carotenoids in studies from chicken through salmon (Ahi, Lecaudey, Ziegelbecker, Steiner, Glabonjat, et al., 2020;Gazda et al., 2020;Lehnert et al., 2019;Toews et al., 2017;Twomey et al., 2020). Over 30 additional candidate pigmentation genes were differentially expressed in red Matadero females relative to dull females of other populations. However, because the Matadero population differed so extensively in expression patterns from the other stream populations, inferring which genes likely play a causal role in coloration is difficult.
In British Columbia's Little Campbell population, csf1ra was significantly elevated in red-throated females, relative to dull females of other populations, but it was not similarly elevated in red Matadero females (although it was in Matadero males). This gene influences carotenoid-based coloration through roles in the development of xanthophores/erythrophores (Kottler et al., 2013;Patterson & Parichy, 2019;Salzburger et al., 2007).

| Tissue-specific and developmental timing of expression of key genes-limitations of current data
A caveat regarding our data is that they are from a single tissue at a single developmental point. Important processes in color pattern development may well take place prior to the stage examined here.
Similarly, in the liver and in other organs and tissues, the expression of key genes may be different from in the skin, and play a role in mediating the relative abundance of the pigments responsible for color patterns (Patterson & Parichy, 2019; but see Gazda et al., 2020). In addition, we did not attempt to evaluate alternative splicing (e.g., Howes et al., 2017).

| CON CLUS IONS
We find that the genes showing correlated expression with female red coloration in stickleback differ extensively between populations. Our results also suggest that the genes mediating variation in female red coloration within populations are largely different from those associated with differences between populations. However, it was unclear from our data whether genes associated with female red chroma were similarly correlated across the sexes, or differed between males and females. Several candidate pigmentation genes from studies of other taxa were identified in our inter-population analyses but very few from analyses within populations. TTC39B, located in a QTL for female red coloration for one of our study populations, showed the most consistent patterns. It will be important to build on our investigation of naturally occurring variation in gene expression with manipulations using CRISPR and related methods, especially for TTC39B.

ACK N OWLED G M ENTS
The

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.